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NOMENCLATURE 

r,  9,  z       cylindrical  coordinates  used  to  describe  the  undeformed 
configuration  of  the  plate. 

h(r),  a  thickness  function  and  radius  of  plate. 

h~  thickness  at  center. 

t  time  variable. 

u»w  radial  and  transverse  displacement  of  the  middle  plane 

er»  eQ  radial  and  circumferential  strains. 

ar>  cr  radial  and  circumferential  stresses. 

all'a12'a22  stress-strain  relation  coefficients. 

Nr>  NQ  middle  plane  forces  per  unit  length. 

Mr>  MQ  bending  moments  per  unit  length. 

Q(S),  Q*(C)  dimensionless  loading  distributions. 

q(r,t)  loading  intensity. 

K  kinetic  energy  of  the  plate. 

Vs'  Vb        strain  energy  due  to  stretching  of  the  middle  plane  and 
due  to  bending  of  the.  plate  respectively. 


W 


c 


work  done  on  the  plate  by  the  external  forces. 


v  Poisson's  ratio  =  -  a19/a99 


ratio  of  elastic  constants  =  a, ,/'a00 

11   22 

D(r)  flextural  rigidity  of  the  plate  =  a22h3/12(a1;La22-a?0) 

ft  <P  stress  functions 

£»  T  dimensionless  space  and  time  variables  respectively. 

X  dimensionless  transverse  displacement 

gCS),  fCO     shape  functions  of  vibration. 

A,  a  amplitude  parameters. 


iv 


X  nondimensional  nonlinear  eigenvalue. 

1/° 

(*)  =  (A.)  '  nondimensional  angular  frequency. 

Y,  Z,  H  (6x1)  vector  functions 

M,  N  coefficient  matrices. 

0  C3xl)  null  matrix 

'r,  't  partial  derivatives  with  respect  to  r  &  t 

n(5)  variable  thickness  function. 

t  frequency  parameter. 

5x  first  variation  of  x  =  SGSincor 

n  adjustable  data  in  the  related  initial  value  problem. 

1  }  indicates,  a  column  vector. 
A  Del  operator. 
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INTRODUCTION 

Composite  materials  find  large  application  in  design  of  structural 
elements  in  the  present  age.   These  structures  which  are  mainly  in  the 
form  of  plates,  are  subjected  to  severe  operational  conditions,  and  should 
thus  be  able  to  withstand  large  amplitudes  of  vibration.   If  the  amplitude 
of  vibration  is  of  the  same  order  of  magnitude  as  the  thickness  of  the  plate, 
then  the  deformation  of  the  mid-plane  can  no  longer  be  neglected.   In  the 
development  of  a  suitable  thin  plate  theory,  anisotropic  properties  and 
geometric  non-linearities  arising  in  the  coupling  of  membrane  and  bending 
theories  should  be  included.  The  resulting  governing  differential  equations 
can  be  solved  by  approximate  numerical  methods,  due  to  the  complexity  of  the 
problem. 

In  1960,  Kazimierz  Borsuk,  determined  a  method  to  solve  in  an  accurate 
manner  the  problem  of  free  vihratlon  of  circular  cylindrically  orthotropic 
plates.   In  1969  A.  P.  Salzman  and  S.  A.  Patel  used  the  method  of  separation 
of  variables  along  with  Frobenius'  method  to  determine  the  frequencies  of 
clamped  or  simply  supported  solid  circular  variable  thickness  orthotropic 
plates.   In  1971,  K.  Vijayakumar  and  C.  V.  Joga  Rao  determined  the  axisym- 
metric  vibration  and  buckling  of  polar  orthotropic  circular  plates.   In  1973 
C.  L.  Huang  and  H.  K.  Woo  used  the  Ritz-Kantorovich  method  to  determine  large 
oscillations  of  orthotropic  annular  plates.   In  1974,  G.  K.  Ramiah  and 
K.  Vijayakumar,  determined  the  vibration  of  polar  orthotropic  annular 
plates. 

The  above  and  many  other  investigators  (18-30.)  have  worked  on  either 
solid,  circular,  variable  thickness,  anisotropic  plates  or  annular  orthotropic 
plates.   This  present  investigation  is  thus  concerned  with,  harmonic,  large 


amplitude,  free  vibrations  of  orthofropic,  axisymmetric,  annular  plates 
of  variable  thickness. 

The  essence  of  the  approximate  method  is  to  approximate  the  continuous 
system  by  a  discrete  one  having  a  finite  number  of  degrees  of  freedom.   The 
discrete  representation  is  achieved  through  an  assumed  space  mode.   Sub- 
stitution of  this  in  the  differential  equations  along  with  the  requirement 
that  some  measure  of  the  error  be  minimized,  the  assumed  space  mods  can  be 
eliminated.   The  problem  thus  reduces  to  a  nonlinear  ordinary  differential 
equation  with  time  as  an  independent  variable.   This  equation  is  similar  to 
a  one-degree  of  freedom  Duf f ings  equation  (5) . 

This  present  work  assumes  the  existence  of  harmonic  vibrations.   The 
time  variable  is  eliminated  by  the  application  of  a  Ritz-Kantorovich 
averaging  method.  The  basic  governing  equations  thus  reduce  to  a  pair  of 
ordinary  differential  equations,  with  a  reformulated  set  of  boundary  condi- 
tions. A  numerical  study  of  these  equations  is  proposed  by  introducing  the 
related  initial  value  problem. 

The  cases  considered  are,  a  parabolic  variable  thickness  annular  plate, 
with  the  variable  thickness  function  of  the  form 

n.  =  Q.  815  -  Q.  5  X2 
and  a  convex  variable  thickness  annular  plate,  with,  the  variable  thickness 
function  of  the  form 

n_=  1.0  -  Q.5  X1'2  . 
Both  the  above  plates  have  the  same  volume  and  the  same  boundary  conditions. 
The  boundary  conditions  are  free  on  the  outside  and  fixed  on  the  inside. 
The  corresponding  curves  for  the  frequency  responses,  bending  stresses, 
and  membrane  stresses  are  presented. 


CHAPTER  I 

DERIVATION  OF  THE  GOVERNING  EQUATIONS 

Consider  a  thin  annular  orthotropic  plate,  the  elastic  properties 
of  which  are  different  in  the  radial  and  circumferential  directions.   The 
fundamental  assumptions  made  as  regards  to  the  flexural  deformations  of 
the  plate  are: 

1.  The  loads  and  deflections  are  symmetric  with  respect  to  the  z 
axis  which  passes  through  the  center  of  the  annulus. 

2.  The  normals  to  the  middle  plane  in  the  undeformed  plate  remain 
straight  and  normal  to  the  middle  plane  in  the  deformed  plate. 

3.  It  follows  the  Hooke's  Law. 

4.  Transverse  shearing  deformations  are.  not  Included. 

5.  The  maximum  thickness  of  the  plate  is  small  in  comparison  to  the 
radius  of  the  plate. 

Keeping  in  mind  the  above  assumptions,  the  following  strain-displace- 
ment relations  are  written: 
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or  in  indicial  notation  as: 


er  *  u'r  +  Tw'r  "  zw'rr  « 


In  the  9  direction: 


£9  m   r   r  3r 


or  in  indicial  notation  as: 
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£e  =  r"-rw'r  «) 


In  view  of  the  orthotropy  of  the  plate  considered,  the  Hooke's  Law 
can  be  written  as: 

(3a) 

(3b) 

From  the  above  two  equations,  it  follows 
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where  a,,,  a„»,  an „  are  the  elastic  constants  and  a  ,  a.  are  the  radial  and 
11   22   12  r   8 

circumferential  stresses.. 


Resubstituting  e.  and  e  from  (1)  and  (2),  we  have 
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Expressions  for  the  radial  and  circumferential  forces  per  unit  length, 
N  and  Nft,  are  obtained  by  integrating  the  respective  stresses  across  the 


thickness  of  the  plate 
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or, 
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where , 


c  =  -   ratio  of  properties  in  radial  and  circumferential 

22   directions. 


12 
v  =  =  poisson's  ratio. 

a22 
The  radial  and  circumferential  moments  per  unit  length,  M  and  Mfi , 


are  obtained  by  integrating  the  moments  of  the  forces  about  the  middle  plane 
across  the  thickness  of  the  plate. 

fh/2  a,,       rh/2 


rn/z 

M  -      a 
r   J -h/2 


zdz 


11 


alla22  -  a12  '-h/2 


"     !l2   ' 


zdz 


(Sa) 


Me  = 


h/2 
-h/2 


a„  zdz  = 


22 


h/2   r     a12 

t- e 


alla22  ~  a12   h/2 


i  9   a22  rj 


zdz 


(8b) 
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The  Energy  Method 

The  extended  Hamilton's  Principle,  which  states  that,  within  the 
interval  of  time,  t,  and  t9 ,  the  first  variation  of  the  action  integral 
is  equal  to  zero,  is  made  use  of  here;  i.e.; 

C2 

<5     Idt  =  0  (10) 


i: 


Here  the  Lagrangian  is  L  =  K-V  -V  +W  (11) 


where 


K  =  Kinetic  Energy 

V  ■  strain  energy  due  to  stretching  of  the  middle  plane 

V  =  strain  energy  due  to  the  bending  of  the  plate. 

W  =  work  done  by  the  time  dependent  external  forces. 
Neglecting  the  radial  part  of  inertia  force,  8u/3t  «  3w/8t,  the 
kinetic  energy  is 
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2.   The  strain  energy  due  to  stretching  of  the  middle  plane  is  obtained 
as  follows: 

I*  N  6    N.e:- 
;c 
Substituting  values  of  N  ,  N  ,  e.  ,  e  and  rearranging, 
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3.   The  strain  energy  due  to  bending  of  the  plate: 
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Substituting  M  and  M„  and  rearranging, 
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4.   The  work  done  by  the  exciting  force  function  p(r,t), 
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Substituting  equations  (12),  (14),  (15)  and  (.16)  into  (11)  we  obtain 
the  Lagrangian  as: 
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Now  the  integral  can  be  written  symbolically  as: 
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by  partial  integration. 

Substituting  values  of 


dr 
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and  t —  into  equation  (21a)  we  obtain: 
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By  similar  method,  substituting  values  of 
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and  - —  in  equation: 
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And  as,  fiv  •  0  it  L ,  t,,  the  first  integral  =  0,  hence  the  above 
equation  reduces  to: 
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And  finally  we  have, 
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after  integrating  by  parts  twice,  and  substituting  the  values  of  the 
various  constituents  of  the  equation  we  have, 
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Substituting  equations  (19)  to  (24)  into  (18)  we  get: 
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For  equation  (25)  to  hold  true,  the.  Integrands  in  the  double  and 
single  integrals  should  vanish  separately.  * 
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The  double  integral  yields  the  Euler-Lagrange  equations: 
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+ +  TW  vt       +  -t— )   +  —  (uw       +  u  w   )  J   +  h   (cu  w 


c     3    .    v  J 


•      o     ~     ™  ■      ~     /      i  \iav»  .      u    „      j  i      i      u     CCU    W       +   —  W       1-    —  UW     J  f 

r         2     r  rr       2r         r         rr         r  r  rrr       2     r       r       r  J 

(26) 


and 


r  cu  w  1 

h-jcu       + ^T  +   (c-v)  -t—  +  cw  w     \ 

{     rr         r  2        v  2r  r  rrj 

j_  u  /  •    c     2    .    vuV       ., 

+  h  -(cu     -r  ■=-  w    +  — >•  =  0 
r[     r       2     r         r/ 

The  single  integrals  yield,  the  boundary  conditions, 


(27) 


w=0   or   -  —  w  (•=■  w  +  cu  +  — )  +  D(r) 
&  r  2  r     r   r 


r     w  cw  -\ 

r  ,  rr  - 

+  cw    + 


2   "  rrr     r 


+  D  (cw   +  —  w  )  =  0  ■+■  shear 
r   rr   r  r 


■  0 


deflection  ■  Q 


f     =  0   or   D(r)(2crw   +  2vw  )  ■  0  +  Moment  =  0 
r  rr     r 


slope  =  0 
and 


cw 


u  =  0   or   rh(r)  (cu  +  —~-  +  — )  ■  0  ->  Force  =  0 


Equation  (26)  can  be  expressed  as: 


13 


D(cw  +  —  w         ttw       +-^-w)+D   (2cw         +   (2c+v)  -  w       -  ^v  w  ) 

rrrr         r     rrr         2     rr         3     r  r         rrr  r     rr         2     r 

r  r  r 

+  D     (cw      +%)  +  ptar       =  p(r,t)  +  |  -  f-  [hrw   (cu     +  |-  w2  +  -  u)  ] 
rr       rr       rr  tt  tsrar  rr       2r       r 

(28) 
Proceeding  further  with  the  stress   formulation,    and  substituting  the 
following  into  equations    (27)   &    (28) 

Y  =  rN  and  |£  =  N 

r  3r  9 

3uJ    1     2       a22    ,M  AT   ,        a22    rV 

37  +  r  wr  =  IT  (Nr  "  VV  =  ~  (7  "  VV 

H  =  !|2    (         _  v 
r         h  r       r 

We  have  the  following  equation: 

D(cwrrrr  +  H  Wrrr  "  \  Wrr  +  \  Wr}  +  DrC2cwrrr  +  (2c+v)  r  Wrr  "  \  Wr 

r       r  r 


+  Drr<cwrr  +  7  wr)  +  Phwtt  =  P(.r,t)  +  |  [*rflr 


(29) 


rr    r     Z     h     r    r    ^a??r  r 


(30) 


and  the  boundary  conditions  become: 


w  =  0   or   D[cwrrr  +  |  wrr  -  i]  +  Dr[cwrr  +  %f]  -  ^f  ^  =  0 


w  =  0   or   cw   +  —  w  =  0 
r  rr   r  r 


¥=0   or   cT-v—   =  (1 
r     r 


Using  the  substitutions, 


w 
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I  =T 


D, 


T  =  t[ 


ha4(l-v2) 


,1/2 


22 

♦  ■  T~  v 


hQa 


D  - 


"1  *3  3 

k    .   r     r  A    ho  n 

12(c-v  )a 


22 


D„  ■ 


0 


0   12a 
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we  get  the  following  non-dimensional  form: 


n3[cX""  +|£x'"  -  ] 

E 


,   X'  *  +  ^  X'  ]  +   (n3)      [2cX' ' '   +   (2c+v)  |  x' ' 


-^2  x']  +  GT>U  [cX"  +|x']  +  n 


r     2i 

c-v    I 


I   x. 


1    2J  A^T 

1-v   ' 


=  I2(c-v2)a99(g-)3  p  +  12(c-v2)    (|-)2  |  [x>]' 

o 


22  vh 


(31) 


and, 


[c*" 


♦^-^♦^W  +  ^+fc 


25   [xC] 


(32) 


The  boundary  conditions  become, 


X1  =  0   or   cX"  +  j  x'  ■  0 


(33a) 


X-  0   or    [ex'"  +  §  x"  -^X'l  +^  [cX"  +^X*3 


-12(.c-v2)  (|-)2  -A-  (x^)  =  0 

o  n.  % 


(33b) 


<{•  *  0   or 


c$'  -  -p-  <J>  =  0 


(33c) 
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CHAPTER  II 


APPROXIMATE  ANALYSIS 


There  is,  at  present,  no  exact  method  known,  for  the  solution  of  the 
differential  equations  C31)  and  (32),  also  the  standard  fourier  analysis 
used  in  linear  vibration  problems  is  not  applicable,  because  the  nonlinear 
character  of  the  differential  equation,  causes  coupling  of  vibration  modes. 

Consequently,  this  nonlinear  coupled  problem,  can  only  be  solved  by 
some  approximate  numerical  method.  Approximate  solutions  of  large  amplitude 
vibrations  can  be  achieved  by  separation  of  variables  method,  or  implementing 
function  space  methods  to  eliminate  the  space  coordinate  with  an  assumed 
mode  shape  function.  The  problem  is  thus  reduced  to  a  non-linear  ordinary 
differential  equations  with  time  t,  as  an  independent  variable.   The  result- 
ing one  degree  of  freedom  Duffings  equation  is  solved  and  the  solutions 
are  in  terms  of  elliptical  functions.   This  is  called  the  assumed  -  space  - 
mode  solution.  The  Kantorovich.  Averaging  method  is  proposed  to  find  an 
assumed  time-mode  solution  of  the  equations  (31)  and  (32)  and  the  boundary 
conditions  equations  (33) . 

Kantorovich  Averaging  Method: 

A  sinusoidal  form  of  the  loading  intensity  is  assumed  here: 

P(£,t)  ■  QCU  Sin  wx  (.34) 

also,  the  steady  state  response  can  be  closely  approximated  by  the  expressions, 

X(£,x)  =  G(£)  Sin  wt  C35a) 

♦  C5,t)  -  F(£)  Sin2o)T  (35b) 

where  G(£)  and  F(5)  are  the  undetermined  shane  functions  of  vibrations. 
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Substituting  equations  (35)  into  equation  (32),  we  have 

■$*tf-£tfc****i-.f*»fr-«         (3«, 

As  the  expressions  (34)  and  (35)  cannot  satisfy  equations  (31)  for 

all  x,  the  integral, 

W*f3[cx?,''+^x,',-7*,,+7*'1 

+  (n3)   [2cx...+i2c+vix„   1  xt] 

+  (n3)^  [cX"  +|-X']  +  n 

r  2. 

c-v 

XTT 

1    2J 

kl-v  J 

-  (c-v2)P  -  12 (c-v2)  (a-)2  |  (X'+)*J  6X£d£              (37) 

where     P  -  12a22  (g)3  p, 

is  used  to  obtain  a  governing  equation  which  closely  resembles  equation 

(31) ,  within  the  limits  of  assumed  form  of  motion  and  loading  as  given  in 

equations  (34)  and  (35). 

• 

Substituting  expressions  (34)  and  (35)  into  (37)  and  equating  the 

average  virtual  work  over  a  period  of  oscillation  to  zero,  or  explicitly: 

f2ir/w 

V     =         I.dT  =  0 

yields: 

n3[cG""  *2|c,.e   1  Qlt  +  1  QI]+  (n3    [2cG,,f  +  C2c+v)  G„ 

r          s                6                   5 

9 

-  K  G'  ]  +  (n\F[cG' '  +  £  G'  ]  -  u2n  ^=M  G 
5                 '       .  Ofc-v 2) 

-  9  (c-v2)  (*-)   |  [G«  fV    =  (c-v2)  Q                     C38) 

o     * 
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The  problem  therefore  becomes  governed  by  a  pair  of  nonlinear  coupled 

ordinary  differential  equations  (36)  and  (.38) .  For  convenience  in  conducting 
a  parametric  study,  let 

G(5)  -  A  g(5)  (39a) 

FCC)  -  A2f  (?)  (39b) 

where  A  is  amplitude  parameter,  and  g(£)  and  f (?)  are  shape  functions. 
Substituting  these  into  equations  (35)  and  (38)  we  get, 

C^*Ca-£5)f  "  t1—^)^*^  (g')2  =  0  (40a) 

en3  g"  "  +  A  «|  +  6cn'n2)g' "  +  C-  4  +  ^^^ 

r  % 

+   3cn,,n2  +  6cn.Cn')2)  .g" 

o 


+  (\  n.  -  3n.'n  ^  +  C3il,,ii  +  6nCn.*)  )£>  g'  - 


nA y~  g 


r  r  *  ci-v ) 


-  9(c-v2)a|  (g'f)  =  (c-v2)  51 


This  can  be  expressed  as: 

Ag""  +  Ag"'  +  Ag"  +Ag*  -  nA-^^-g  -i^Lla(g.  f). 

(c-vZ)  ^ 
=  — ~   Q*  (40b) 

where 

3 


A1  =  en" 


5 


A„  =  —  n.  +  6cnTrL 


2 
A3  =  -  ~  a3  +  3C2c-h;)  3^Jk-  +  3cn"n.2  +  6c  Cn.')2  n 
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2 
*    1  „3   3n,'n    /,  ti  2   ,  ,  ,,2.  v 

c  r  * 

a  =  A(-) 
o 

.  _  2 

A  =  CJ 
O 

The  above  equations  together  with,  the  boundary  conditions  selected 
from  Table  I,  constitute  a  two-point  boundary  problem  which  is  solved 
through  the  solution  of  the  related  initial  value  problem. 


Table-I 
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Type  of  Edge 


Boundary  Condition  at  Edge 
£,  =  R  or  1 


Clamped 
Immovable 


g  -  0 
g'  =  0 


cf'-if-fl 


Clamped 
Movable 


g  =  0 

g'  .  o 


f-0 


Hinged 
Immovable 


g  ■  I 
eg 


"  +  -  - 


g'  =  0 


cf»  -i  f  -0 


Hinged 
Movable 


g  =  0 

eg"  +7-g'  =  0 


f- 


Free 


eg"  +7g'  "  0- 


eg' 


-[V^flg'-O 


f=o 
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CHAPTER  III 


NUMERICAL  ANALYSIS 


The  solutions  of  nonlinear  boundary  values  and  nonlinear  eigenvalue 
problems,  are  very  complicated  and  hence  these  are  solved  by  converting 
them  to  initial  value  problems. 


Initial  Value  Method 

Due  to  the  extremely  nonlinear  form  of  these  equations,  after  conver- 
sion to  an  initial  value  problem  the  shooting  technique  is  used.   An 
associated  variational  problem  is  developed  and  used  in  Newton-Raphson 
iteration  scheme. 

The  governing  equations  (40a)  and  (40b)  can  be  written  as  a  system 
of  six  first  order  differential  equations, 


H  =  H(C,  Y,  a,  A,  Q*) 


R  <  £  <  1 


(41) 


where 


g 

yl 

g* 

?2 

g" 

^3 

Y(£)   =  '- 

«.  »     1     f 

=5   . 

o 

T4 

f 

*> 

V 

7*J 

C  )'  =~   ,  and 


H  is  the  appropriately  defined  (6x1)  vector  function: 
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H  =  ' 


A, 


nx 


A.,    y4 


y /.  "  a    y^  ~  a    y^  "  a    y 


(c-v2) 


V 3 


A      *2        A      Jl  2 


9 (c-v2)      ,  .  .         (c-v2)   _ . 

^A~  a(^3y5  +  W6>  +  7T~  Q 
1  /a  A.. 


y6 


,,        n.'      _.    y5        ,_        n*    _.    y6  n,      ..2 

(1 .  _  v5)  _ .  a  _  _  0  _  _  _  <y2) 

The  parameters  a  and  X  are  at  present  not  known  and  hence  two  additional 
restraints  are  imposed  to  evaluate  these.   One  component  of  Y(l)  is  normalized 
to  fulfill  the  requirement. 

The  boundary  conditions  can  be  expressed  as: 

w 

MY(i)  -  i° 
Q 

PS 


(42a) 


NYCR)  =  - 


0 


(42b) 


where,  M&N  are  (4x6)  and  (3x6)  coefficient  matrices  of  rank  4  and  3 
respectively.   The  first  row  of  M  normalizes  a  component  of  Y(l),  and  the 
remaining  rows  of  M&N  are  obtained  by  taking  into  consideration  the  boundary 
conditions  at  the  two  ends. 

The  corresponding  initial  value  problem  may  be  expressed  as 


dZ 


~  =  HCS,  Z;  a,  \,   0*) 


(43a) 
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Z(l) 


(43b) 


where     y   is  a  (6x1)  vector  of  initial  values. 

Substitution  of  these  initial  values  Z(l)  =  y  into  the  equation  (^2a) 
yields: 

0 


MZ(1)  =  My  =  ■ 


(44) 


As  M  is  of  rank  four,  two  additional  values  are  required,  and  by  the 
implicit  function  theorem, 

Y  ■  Y*^,^) 

is  a  solution  of  equation  (44)  ,  and  r]   ,   r\     are  arbitrary  initial  values. 
The  related  initial  value  problem  can  thus  be  written  as: 


||=  au,  z;  a,  x,  q*) 


_* 


Z(i)  =  Y  (n.vr)2l 


C45a) 


(45b) 


The  above  contains  the  initial  values,  which  satisfy  the  boundary 
conditions  (42a) . 

Assuming  a  continuous  function  Q*(£),  a  solution  of  the  initial  value 
problem  (45)  is  obtained  over  a  closed  interval  [R,l],  and  is  denoted  by 
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r 

-     -      -  -    P1 

Z  =  Z(l;  n,  a)    ,   n  =jn2  • 

[*  j 

From  the  boundary  condition  (42b), 

NZ(R;  n,  a)  ■  0  (46) 

Stating  a  well-know  matrix  theorem,  "For  a  system  of  equations 
NZ(£;  n,  a)  =  0  a  necessary  and  sufficient  condition  for  a  unique  solution 
n  "  n(a),  is  that  the  determinant  of  the  Jacobian  matrix,  J  =  —  [NZ(R,  it,  a)] 
is  not  equal  to  zero,"  assuming  also  that  Z  is  continuously  differentiable 
with  respect  to  n  and  a. 

Thus  there  exists  a  locally  unique  function  at  £  =  R,  such  that, 

NZ(R,  n.(cO,  a)  =  0 
Or,  we  can  put  it  as 

Y(£,ct)  =  Z(g,  a(a),  a)  . 

This  forms,  an  a-dependent  family  of  solutions  to  (42)  each  of  which 
is  a  solution  to  the  initial  value  problem. 

For  a  fixed  value  of  a,  say  <x°,  equations  (46)  reduce  to  three 
trancedental  equations, 

NZCR,  5,  a°)  =  0  (47) 

A  root  n.  of  (47)  may  be  obtained  by  Newtons  iteration  method.   Starting 
with  an  initial  guess,  n.  3  rL  the  iterative  sequence, 

Vl  =  nk.+  A\   '   k  =  1»2»3»—  (48a) 

is  generated. 

This  can  be  expanded  in  the  Taylor's.  Series.  Retaining  only  the 
linear  terms,  gives, 


.. 
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Ank  -  -[N  ^-  Z(R,r,k,  a0)]"1  NZ(R;  nk>  a°) 


(48b) 


3nk 


th 


where,  at  the  k   step,  the  (6x3)  matrix  J.,  is  defined  as, 


(v  -  (41    - 

1      We-R 


3Z 


PL) 


t'V** 


i  a  1,2,. ..,6 
J  -  1,2,3 


(49) 


to  n. 


Physically,  this  represents  the  change  of  final  values  with  respect 

The  expression  NZ(£,  n,  ,  a  )  represents  the  k   error  vector. 

If  the  initial  guess  n.,  is  in  the  neighborhood  of  n_  ,  then  the  con- 
vergence of  the  sequence  n.  to  the  root  n.  is  feasible. 

In  order  to  generate  the  sequence  n,  ,  it  is  necessary  to  evaluate  the 
matrix  (J*)k  at  each  step,  k,  of  the  iteration  process.  To  do  this,  an 
associated  variational  problem  is  introduced. 

Formally  differentiating  (45)  with  respect  to  n.,  gives 


dC   [3nJ 


3Z 


an 


3K 


W 


3Z 
3ti 


(50a) 


'all  ay* 

3nJc=1  "  3n. 


(50b) 


which  constitute  eighteen  first  order  equations.,  and  a  corresponding  set  of 
initial  values. 

d£     2 

dZ„ 
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dZ3 
dC~=  Z4 


dZ4    2  2  (c-v2) 

dT  =  n  w  ~^2^  yl  "  A4y2  "  A3y3  "  A2y4 


2N   a 


+  9Cc-v  )  |  (7373  +  y2y6)  +   (g^2)   Q» 
Al  i^A, 


(51) 


dZ5 
dT"=  z6 


dZ( 
dT 


.J-n-fwVyj-tt.f.Cl*^- 


2£c  y2 


Differentiating  the  above  wxth  respect  to  (ft.  ,  o_,  X)  we  get  the 
following  variational  equations: 


d 

dC 

f3Z  • 

3Z2 

d 

d£ 

f9V 

SI 

3Z3 

3a, 

d 
dC 

r3Z3> 

_ 

3Z4 
3a, 

d 
d£ 

r3Z.> 

4 

- 

f  2  2 
ft  u 

,        2.    3Z.               3Z„ 
(c-v   )        1                  2 

(1-v2)    3nl          4  3nl 

+  9Cc- 

2v    a 
-v)? 

r3Z                 3Z 

7      4- 

[3H,     5        3a, 

3Z? 

w 

W  A! 

d      1 

-3Z-j 

3Z, 
6 

d£    1 

3tiJ 

an. 

» 

3Z,      3Z. 

A3  3a    2  an, 


3Z, 


Z3  +  in7  Z2 


(52a) 
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2 

dC    I3TL 


i  3Z, 

n         3n, 


/  5 


c£     2 


3Z, 

3n^ 


d_ 

dC 


d_ 
dC 


d_ 
d£ 


f9Z, 


I3  ^2 
f3Z2) 


[3n2J 

3Z 


3n~j 


3n„ 


5 

3ru 


3Z, 


dC  (an2J 


2  2    (c-v2)    9Z1  3Z2  9Z3  3Z4 

V35T    A4^"A3^"A2^ 


CL-v   )    °"2 


+  9(.c-v2)  | 


3Z  3Z 

Z,  +  t-5-  Z.  + 


3^ 


3tW 


3Z2 
"Sn^  Z6 


3n2  S 


]}.". 


(52b) 


dS     3n0 


3Z, 


an. 


d_ 
d5 


,  3ZC 

1  Sr^ 


3Z, 


/c5  -[1-r5]^' 


3Z, 


£c     2   3r 
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■3Z. 


3X 
3Z. 


d_ 

d_ 

dC 


d£  [3X 

r3Z4 


3X 
3Z, 


3Z, 
'§T 

3Z, 

w 

3X 


d_ 
dC 


3X 


2^  3Z       3Z       3Z.      3Z. 

_i  _  A  £  _  A  1  _    4 

c   2}  3X    A4  3X    A3  3X    A2  3X 


2  2  Cc-v  ) 


+  nTZ,  +  9Cc-v~)  =■ 


3Z. 


C  |3X 


A.  + 


92, 


C52c) 


3Z 
+  3X 

2                9Zfi 
LZ6+3X6Z 

d 

r3V 

az, 

Q 

dC 

[ax  ] 

ax 

Si 

r3Z6l 
lx  J 

-  El  - 

31  *»,  "5 

n    *  J  3X 

_  il_ 

3Z 

Z„    7T^ 

=1)  '  Al 


9  «    9Z< 

/^-^-^ar 


/ 1 


Cc  2  ax 


For  a  given  vector  ri_  and  a  =  a°,  this  derived  problem  along  with,  the 
initial  value  problem  (.45)  may  be  integrated  simultaneously  on  the  interval 
[R,l].   Corresponding  to  a  given  value  of  n.  and  a  =  a°,  the  calculation  of 
the  resulting  solution  to  the  variational  problem  at  £  =  1  provides  the 
Jacobian  (J-J.   By  setting  n  =  a,  and  integrating  equations  (45)  &  (50) 
from  C=l  to  £=R,  gives  the  first  correction  vector  rL.   By  repeating  this 
procedure,  the  desired  sequence  rfk  is  obtained,  which  converges  to  n°  within 
a  specified  error  bound  to  the  accuracy  of  the  system. 

Having  obtained  n,  ,  corresponding  to  a°,  the  value  of  a  can  now  be 
perturbed, 


a  =  a  +  Aa°  ■  a 
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The  problem  is  reinstated,  for  this  value  of  a,  starting  from  n_  ■  n  . 
If  Aa  is  small,  then  n.  is  contained  in  the  new  contraction  domain  of 
Newton's  method,  the  iterations  converging  to  n  corresponding  to 
a  =  a  .   Successive  completion  of  this  operation  j  number  of  times,  yields, 

n"  =  n.   Ca  )    ,   i  ■  Q,l,2,...  ,j 
By  setting  a  =  a  +  Acr  =  a    and  starting  integration  from  n  -1 , 

one  obtains  n     provided  Acr  results  in  convergence. 

The  range  of  a  is  limited  as  the  elastic  plate  cannot  withstand 

unbounded  amplitudes. 
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CHAPTER  IV 


NUMERICAL  COMPUTATIONS 

The  above  theoretical  analysis  suggests  the  use  of  a  numerical  inte- 
gration technique. 

Use  is  made  of  a  fourth  order  Runge-Kutta-Gill  method  to  integrate 
the  initial  value  problems  (45)  and  (50)  over  the  interval  [R,l].   The 
following  approach  is  suggested. 

The  problem  is  first  reduced  to  that  of  a  free  vibration  by  setting 
Q*  =  0  and  a  =0.   By  subjecting  this  equation  to  a  particular  set  of 
boundary  conditions  the  linear  eigenvalues,  and  mode  shape  functions  are 
determined. 

This  information  leads  to  a  basis  for  making  a  reasonable  starting  guess, 
n-p  required  by  the  initial  value  method. 

For  n.  =  n-p  the  initial  value  problems  (45).  and  (5Q)  are  integrated 
numerically  over  [R,l]  with  a  step  size.  Au  ■  1/40.   Successive  correction 
is  carried  out  till  all  equations  in  (471  satisfy  the  range  of  prescribed 
error;  this,  being  consistent  with,  the  order  0  (.|Au|4).,  of  Runge  Kutta  Gill 
method . 

By  gradually  incrementing  the  value  of  a   and  restarting  the  correction 
and  integration  procedure  from  the  values  of  (ju,  ru,    X) ,  obtained  from  the 
solution  corresponding  to  the  previous  a,  the  resonance  curve  and  other 
solutions  are  evaluated.  This  procedure  is  terminated  at  a  particular 
value  of  a  ,  because  of  reasons  mentioned  earlier. 
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Cases  Considered: 

The  cases  considered  are: 

(1)  Annular  circular  plate  with,  convex  variable  thickness  and  free 

on  the  outside,  fixed  on  the  inside. 
(.2)  Annular  circular  plate  with  parabolic  variable  thickness  and 

free  on  the  outside,  fixed  on  the  inside. 
The  figures  pertaining  to  the  above  two  cases  are  as  shown  in 
Appendix  I. 

The  governing  equations  and  boundary  conditions  are  written  as: 
dY 


where 


d5 


=  HCS,  Y;  a,  X,  Q*) 
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The  related  initial  value  problem,  is  defined  as, 

dZ 


dC 


=  HCC.  Z;  a,  X,   Q*) 


(54a) 


Z(l)  = 
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and  the  variational  problem  is 
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The  unification  of  the  above  is  as  symbolized  in  equation  (50)  with 
i  ■  (n,,  0^,  X)  while  the  value  of  a  is  held  constant. 


*1 


\r\2 


An. 


3Z 


1 

3a, 

3T12 

3Z2 

3Z2 

3\ 

3n2 

?Z6   v 
3^   5 

3Z5 

3Z6   v  8Z5 
C  3n2   c  3n2 

3Z1 
3X 

-1  *T 

3Z2 

3X~ 

< 

3Z6 

3X 

3ZC 

v   5 

5=R  I 

CZ6  "  |  Z5 


jS=R 


provides  the  linear  correction  of  the  estimated  values  Ol,  ,n9,X),  where 


Z  are  components  of  Z. 


For  each  value  of  a  ,  a  sequence  which  defines  discrete  values  of  y, 
successive  corrections  of  (n1,n2,X)  were  performed,  till  the  final  values 
of  Z(R)  are  satisfied, 


max 
Ki<3 


l^    Z.  «> 


<  0.1  x  10 


-5 


(57) 


where  n. .  =  N. 

Perturbing  the  amplitude  a,  the  process  is  started  using  the  values  of 
ri  obtained  after  the  first  cycle  is  completed.   At  least  five  to  six 
iterations  were  required  for  most  values  of  a. 
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Stresses: 

From  all  the  discussion  carried  out  so  far,  it  is  obvious  that  the 
amplitude  influences  the  distribution  of  bending  stress  to  a  greater  extent, 
as  these  are  related  to  the  derivatives  of  the  transverse  shape  function 

1(E). 

The  expressions  for  the  bending  and  membrane  stresses  are: 
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These  are  the  radial  bending  stress,  circumferential  bending  stress, 
radial  membrane  stress  and  circumferential  membrane  stress  respectively, 
in  terms  of  the  dimenslonless  deflection,  x»  and  stress  function  f, 
respectively.   Taking  the  previous  assumption  into  consideration,  i.e., 

xOS.tJ  =  AgC£)    Sin  cot 

♦  CS,t)   =  A2.f CO    Sin2o)T 

and  also  taking  into  consideration  the  fact  that  when  time,  t,  is  equal  to 

the  odd  multiple  of  ir/2w,  we  have  the  maximum  stresses, 

b  2 
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The  above  equations  were  used  in  the  computer  program. 
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CHAPTER  5:   CONCLUSIONS 

This  study  is  based  on  the  supposition  of  harmonic  oscillations.   The 
assumed  solutions  (35)  contradict  the  inseparability  of  modes  in  Von 
Karman's  dynamic  equations.  Nevertheless,  for  moderate  amplitude  of  vibra- 
tions, physical  arguments  may  be  made  to  justify  such  assumptions. 

The  time  coordinate  function  is  assumed  and  eliminated  by  a  time 
averaging  method.   By  elimination  of  the  time  variable,  an  infinite  number 
of  degrees  of  freedom  in  the  space  coordinate  function  is  achieved.   By 
the  numerical  integration  technique  used  the  solution  of  the  continuous 
system  is  obtained  at  a  number  of  discrete  points.   This  reduces  the 
number  of  degrees  of  freedom  to  the  number  of  points  considered. 

The  two  cases  studied  are  an  annular  plate  with  parabolic  variable 
thickness  and  of  convex  variable  thickness.   Both  are  of  the  same  volume 
and  have  the  same  boundary  conditions  of  free  on  the  outside  and  fixed  on 
the  inside. 

The  responses  of  the  plates  exhibit  a  behavior  similar  to  that  of  a 
hard  spring. 

The  parabolic  variable  thickness  plate  is  stiffer  than  the  convex 
variahle  thickness  one,  as  is  evident  from  the  frequency  responses  obtained. 

The  bending  stresses,  of  the  first  plate  are  slightly  higher  than 
those  of  the  second  plate  and  the  membrane  stresses  are  just  the  reverse. 

The  membrane  stresses  have  significant  magnitudes  even  at  relatively 
low  amplitudes.  This  is  due  to  a  stress  concentration  factor  at  the  edge 
of  the  hole,  and  is  called  the  boundary  layer. 

The  results  obtained  were  compatible  with  those  obtained  by  Sandman 
[2].   In  this  study  the  higher  modea  and  stability  of  vibration  have  not 
been  considered.   So  also,  various  other  boundary  conditions  are  possible, 
these  are  thus  left  open  for  future  investigation. 
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Appendix  A 


Figures 
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Fig.  1.   Circular  Plate  and  the  Polar  Coordinate  System. 
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Fig.  3.   Normalized  frequency  responses 
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Fig.  4.   Shape  function  for  annular,  convex-variable 
thickness  plate. 
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Fig.  5.   Shape  function  for  annular,  parabolic,  variable 
thickness  plate. 
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Fig.  6.  Radial  bending  stresses 
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Fig.  7.   Circumferential  bending  stresses. 
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Fig.  9.   Circumferential  membrane  stresses. 
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APPENDIX  B 


Computer  program  for  annular  orthotropic 

convex  variable  thickness  plate  — 

Backward  shooting 
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$JQB 

C***********# #******************************************* 

C      INITIAL  VALUE  METHGO  -  FREE  VISRATICN  OF  AN  ANNULAR 

C      ORTHOTRCPIC  , CONVEX  VARIABLE  THICKNESS  PLATE,  WITH 

C      BOUNDARY  CONDITIONS  AS  FIXED  CN  THE  INSIDE  AND 

C      FREE  CN  THE  OUTS  ICE 

C************ ****** ************************************** 

C************ ********* ******************* **************** 

C      S=RATIC  CF  ELASTIC  CONSTANTS 

C      E=PCISSCNS  RATIO 

C  QL=UNIFCRM    LCAOING    INTENSITY 

C  A=AMPLITUDE 

C  R=RATIO    CF    INNER    TO    CUTER   RAOIUS 

C      P=EIGE.\VALUE 

C      H*STEP-  SIZE 

C      ETA=TH1CKNESS  FUNCTION 

C      DETA=FIRST  DERIVATIVE  OF  ETA 

C      DDET=SECCND  DERIVATIVE  QF  ETA 

C******* *********** ****************** ******************** 

1  IMPLICIT  REAL*8(/-H,C-Z), INTEGER(I-N) 

2  DIMENSION  ETA (41 ) , XX( 41 ) , Y( 24) ,C( 24) ,TP(3 ,3) , 0 (6, 41 ) 

3  DIMENSION  C  (3  )  ,LVi  (3  )  ,MW  (3  )  ,  ER  (3  ) 

4  DIMENSION  RBS (45)  ,CBS (45) ,RMS ( 4 5) ,CfS (45) 

5  112  F0RMATt5X,'AMF=',C22.14,3X,'FREQ=«,D22.14,3X,«FRER=« 

ID22.14) 

6  113  F0RMA?'9X,'W ,19>, »DW • , 18X, »DDW • , 1 7X , '00CW« ) 

7  114    F0RMAT(4C22.14) 

8  115    F0RMAT(//9X,»F» .lSXt'DFM 

9  117    FORMAT (1H    ) 

10  120  F0P.MAT(5X,«S=«,F10.3,5X,»E=«,F10.3,5X,'QL='  ,F10.3) 

11  121  F0RMAT(//9Xt'STA',l9Xt,PRCF») 

12  122  F0RMAT(6X,'R8S',iaX, 'C3S' , 18X, • RMS • , 18X, 'CMS') 

13  123  F0RMAT(5X,« ITER=» ,12) 

14  S=0.5 

15  H-1./40. 

16  LL=41 

17  KK=9 

18  JK=LLU-KK 
IS  IK=1 

20  CL=O.C 

21  A=0.0 

22  E=l./3. 

23  R=0.2C-0 

24  DA=0.25 

25  P*8.3**2 

26  ETA1=1.4482 

27  ETA2^-0.2426 

28  BE=0.5 

C******* ********************************* **************** 

C       CCNSTRLCT  INITIAL  VALUES 

C********** ************************** ******************** 

29  5CC  IT=l 

3C  00  9  1=1,24 

31  9  Y(I)=C.OC-0 

32  Y(l)=*1.00-0 

33  V(2)=ETA1 

24  Y(3)*-(E*Y(2))/S 

35  Y(4)=((1.*E)*Y(2))/S 

36  Y<6)=ETA2 

37  Y(8)=1.0C-0 


51 


41 

600 

42 

43 

623 

44 

45 

620 

4* 

47 

48 

49 

50 

51 

615 

£2 

624 

38  Y<9)«- E/S 

39  Y(10)-(1.0*E)/S 

40  Y118J-1. 

C*#******** ******** ************************************** 

C       X=INCEPEN0ENT  VARIABLE 

C      INTEGRATION  BY  BACKWARD  SHQCTING 

C ************ **************************** **************** 

X=1.00-0 
00  623  1*1,24 
Q(IJ=O.CC-0 
DO  620  1=1, ft 
DC  I,LL)=Y< I) 
HR=-H 

DO  624  J=2,JK 
K=LL+1-J 

CALL  RKGU,HR,Y,C,P,£,S,E,QL, BE) 
DO  615  L=l,6 
0(L,y)=Y(L) 
CONTINUE 
C************* ********* ************** ** ** **** ************ 

C      ER(I)=ERPOR  VECTCR  FCR  BOUNDARY  CCNCITICNS 
C******** ******************************** **** ************ 

53  £RC1)=C(1,KK) 

J4  ER(2)=0<2,KK) 

55  ERC3)=S*0(6,KK)-<E*0<5,KK))/R 

C ******************************************** ************ 

C      CCNSTKLCT  ERRCR  NORP 

C************ **************** ** ****** **** **************** 

56  DO  26  1=1,4 

57  DER=DABS(ER(IJ) 

58  IF<0ER.GT. 0.10-05)    GG   TO    28 

59  26   CCNTINLE 

60  GO   TC    SCO 

61  28    CONTINUE 

£****** x*4it ************************************** ******** 

C  NEhTONS    METHOO    (ERRCR    NORM=(B)Y=0 

C  TP(I,J)=   Tt-E    JACCBIAN   OF   THE    INITIAL   VALUES 

C        ■     C(I)=CCFRECTICN    VECTCR 

C**** ******************************************** ******** 

62  TP(l,l)=Y(7) 

63  7P(2,1J=Y<8> 

64  TP(3,1)*S*YU2)-(E*Y(11))/R 

65  TPU,2)*YU3) 

66  TP<2,2I=Y<14) 

67  TP(3t21*S*Y(18)-(E*YC171)/R 

68  TP<1,3)*Y<19) 

69  TP<2,3)=Y(2C) 

70  TP<3,3)=S*Y(24)-(E*Y(23)J/R 

71  DET=O.OD-0 

72  CALL  CMINV(TP,3,CET,LW,MW) 

73  00  75  1=1,3 

74  C(I)=0.0 

75  DO  75  J=l,3 

76  75  C(I)=C(I)-TP(I,J)*ERU) 

C************* ******** *********************************** 

C      CORRECTEC  VALUES 

C ********** ********************************** ************ 

77  DO  76  1*1,6 
76       76  Y(I)=0(I,LLJ 
79          Y(2)=Y<2)*CU) 


52 


tO  Y(6)»Y(6)*C{2) 

61  P=P*C(3> 

82  Y«3J  —  (E*Y(2))/S 

e3  Y(4)«((1.+E)*Y(2))/S 

64  00  77  1=7,24 

65  77  YU)=0.00-0 

66  Y<81=1.0D-0 

67  Y<9)=-E/S 

68  Y(i0>=(l.0«-El/S 

69  Y(18)=1.0 

90  IT=IT*1 

91  IFU7.GT.10l  GO  TO  550 

92  GO  TO  tCC 

£*********** *********** ************** ** ****************** 

C  FINAL   RESULTS 

C************ ****** ****************** ************ ******** 

93  9C0    SRA=CSORT(A) 

94  SP=DSQRT(P) 

95  00   795    J=KK,LL,4 

96  DJ=J-l 

97  XX(J)=OJ*H 

98  ETA< J)«1.-BE*<XX(J)**C0.5)) 

99  IF(XXU).GT.O.O)    GO   TO   905 

100  R8SUJ  =  SRA*C(3,J)/2.*<1.-E) 

101  CBS(J)=RBSU) 

102  RKS(J»=A»0(6,J) 

103  CMS(J1=R*S(J) 

104  GC    TO   795 

105  905   RBS<J)=SRA*ETAU)*{S*C(3,J)  +  E*0<2,J)/XX(J))/2.*(S-E**2) 
1C6  CBSU>=SRA*ETA<  J)*(C(2,J)/XX<JME*Dt3,J) )/2.*( S-E**2> 
107  RMS(J)-A*0(5,J)/ETA(j)*XX(J> 

1C8  CMS(J1=A*0(6,J)/ETA( J) 

1C9  795    CCNTIME 

£****** ******************************  ******************** 

C  FOP  FREQUENCY  RATIO 

C ************ ******************************** ************ 

110  lF(A.GT.O.O)  GO  TO  906 

111  SPC=SP 

112  9C6  SPR=SP/SPO 

113  WRITE<6,117) 

114  WRITE(6,120)  S.E.CL 

115  WRITE<6,117J 

116  WRITE(6,112>    SPA.SP.SPR 

117  WRITE16»U71 

118  WRITE(6,113) 

119  DO    901    J=KK,LL,4 

120  901    WRITE<6,114>    (0(  I, J), 1-1, 4} 

121  WRITE(6,115I 

122  DC    902    J=KK,LLf4 

123  902    WRITE(6,114J     (0(L,J)  ,L*5 ,6) 

124  WPITE(t,117) 

125  WRITE16.122) 

126  00  902  J=KK,LL,4 

127  903    WRITE<6,114)    RBS< J) ,CBS< J S ,RMS< J> tCPS( J) 

128  WRITE(6,117J 

129  MUTE(6,123l    IT 

130  WRITE(6,121) 

131  00    924    J=KK,LL,4 

132  924   WRITE(6,114»    XX< J) ,ETA( Jl 
123  WRITE(6,U7) 
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c*«,*»****************************** ********************** 

r  PEBTUKtATIGN  CF  AMPLITUCE                    ...—** 

c******************************************************** 

134  A*A+CA 

125  IMIK+1 

l->6  IFUK.GT.26I  60  TC  550 

137  ETA1=C(2,LL) 

136  ETA2=C(6,LL) 

12$  P»(SP-0.3)**2 

140  GO  TO  5C0 

141  550  STCP 

142  END 

143  SUBROUTINE  PKGtX ,H,Y ,Q,P, AP,S , E,QL, BE) 
,44  IMPLICIT  REAL*8<A-H,0-£JtINTEGEP(I-N> 

145  0IMENS10N  Y(2A)tC(24),0Y(24),A(2) 

146  A(1)=C. 2928932188135 

147  A<2J=1.1071C678lie65 

148  H2*0.5*h 

149  CALL  DERIV(X,H,Y,DY,P,AP,S,E,CL,BE) 

150  00  13  1=1. 24 

151  R  =  H2*CYU)-CII) 

152  Y(I)=Y(I)+R 

153  13  Q(l}=C(l)+-3.0*R-H2*OY(I) 

154  X=X*H2 

155  00  60  J  =  l,2 

156  CALL  CERlVlX,H,Y,CY,P,AP,S,EfQL,BE) 

157  OC  20  1=1,24 

158  R=A(J1*(H*0Y(I)-C( !) J 

159  Y(IJ=Y(I)+R 

160  20  Q(I)=Q<1)«-3.0*R-A(J)*H*DY(I> 

161  60  CONTINUE 

143  CALL  OERIV(X,H,YTOY,P,APtSfEfCLfBE) 

164  DO  26  1=1,24 

165  R=<H*DYCIJ-2. 0*0(1)1/6.0 

166  Y(H  =  Y(I)+R 

167  26  Qm=Q(I)-»-3.0*R-H2*CY(H 

168  RETURN 

169  ENO 

170  SUBROUTINE  DERI  V  (X  ,h  ,Y  ,0Y  ,  P,  AP,  S,  E,  CL,  BE) 
l7i  IMPLICIT  REAL*8(A-h,C-I),INT£GER(I-N) 

172  OIMENSICN  Y124),CY(24) 

173  ETA=1.-BE*(X**<0.5)  ) 

174  D£TA=-0.5*BE*(X**(-0.5)) 

175  D0ET  =  0.25w8E*(X*'»(-l-5J) 

176  DO  10  1  =  1,3 

177  10  DY«I)=Y(I*1J 
176  DY(5)*Y(6) 

179  DO  12  1=7,9 

180  12  0YU)=Y(H-l) 

181  DY(11)=YC12) 
162  00  15  1=13,15 

183  15  DYCIJ-YU  +  l) 

184  DY(17)=Y(18) 

185  00  16  1*19,21 

166  16  DY(I)*Y(I+1) 

167  DY<23)=Y<24) 

188  IFIX.GE.0.1C-02J  GO  TO  17 
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18*  0Y(4)*3.»P*Y{l)/8.*<ETA**2}*27.Ml.-E**2)*AP*Y<3)*Y(6} 

1/4.*<ETA**31 
1*0  0Y<4)=0Y(4)-9.*(1.+E)*0CET+Y(3)/8.*ETA 

1*1  IF(AP)     18,18,19 

152  19    DY(4)=0Y(4)+3.»(l.-E**2)*Ql./8.*{ETA**3}>>OSaRT<AP) 

153  18   DY(4)=0Y<4> 

1**  DY<10>=3.*P*Y<7>/8.*(ETA**21>27.*(1.-E**2)*AP*{Y(9)* 

1Y(6)*Y(3)*Y(12)]/4.>META**3I 
1*5  DY(10)=0Y(10)-9.*(1.+E)*00ET*Y(?)/8.*6TA 

196  DYU6)=3.*P*Y<I3)/8.*CETA**2)*27.*(  l.-E**2)*AP* t Y( 15)* 

1Y(6)*Y(18)*Y(3) )/4.*(ETA**3) 
1*7  0Y(16)=CY(16J-9.*(1.4E)*0D5T*Y( 15)/6.*ETA 

1*8  DY(22)=3.*P*Y(19)/8.*JETA**21+3.*Y(U/8.*(ETA**3)+27. 

1*<1.-E**2)*AP*(Y<21)<Y(6)+Y(3)*Y<24>)/4.*<ETA**3) 
1*9  DY(22J=0Y<221-9.*{1.+E1*0DET*Y(21)/8.*ETA 

200  0Y(6)=0.0 

2C1  DY<12)=0.0 

202  DY(18)=0.0 

203  DY<24)=0.0 

204  GO    TO    70 

205  17   B0*(9.*(S-(E**2) >  >/ (S*X*(  ETA**3  ) ) 
2Ct  Bl=(S-E**2)/< ll.-E**2)*(S*(ETA**2))  ) 
207  B2=((6.*CETA)/ETA)+(2./X) 

2C8  B3  =  (-1./(S*(X**2)))+(3.*DET/!*<2.*S  +  E)  )/ <S*X*ETA)+(3 .* 

1D0ET)/ETA-K£.*(CETA**2)1/(ETA**2) 
2C9  B4=l./<S*(X**il>-(3.*CETA)/<S*tX**2)*ETA)f(3.*E*DDE7>/ 

1(ETA*S*X1*6.*E*(DE7A**2 ) / ( S*X* ( ET A**2 ) ) 

210  DY<4)*BC*AP*( Y<3)*Y{5)+Y(2)*Y(6)>+B1*P*Y(1)-84*Y<2) 
1-B3»Y(3J-B2+Y(4) 

211  1F(AP)     100,100,101 

212  101    0Y(4)=CY(4)+(S-(E**211*QL/(S*tE7A**3)*0SaRT(AP») 

213  100   DY(4)=0Y(4) 

21*  DY<6)=(1.-E*X*DETA/ETA)*Y(5)/(S*<X**2))-(1.-X*DETA/ETA 

1J*Y16)/X-(ETA*(Y12)**2))/(2.*X*S) 
215  DY( 10)=8C*AP*<Y(fc)*Y<8H-Y(5)*Y(9>+Y<2>*Y(121+Y(3)*Y(ll 

1)  )*81*P*Y(7)-B2*Y( 10)-83*Y(9)-B4*Y( 6) 
21<  DY(12)=(1.-E*X*0ETA/ETA)*Y(U )/ IS* (X**2 ) )-{ l.-X*OETA/ 

1ETA)*Y( 12)/X-ETA*Y<2)*Y(3)/(S*X) 

217  DY(16)  =  eO*AP*(Y(5)*Y(15)+Y(3)*YU71+Y(  14  J *Y (6) +Y< 2)* 
1Y(18) )+ei*PfY(13  1-B2*Y(16)-e3*Y(15)-B4*Y(14) 

218  DY(181=<1.-E*X*CETA/ETA)*Y<17)/<S*(X**2))-<1.-X*DETA/ 
1ETA)*Y(18)/X-ETA*Y(2)*YU4)/(S*X1 

21*  0Y(22)=B0*AP*(Y(21J*Y(5)*Y(3)*Y(231*Y(20)*Y(6)4Y(2J* 

1Y(24J )*ei*P*Ytl9]*Bl*Y(l)-B2*Y(22)-B3*Y(21)-B4*Y(20) 

22C  DY(24)=(1.-E*X*0E7A/ETA1*Y(23)/(S*(X**2))-(1.-X*DETA/ 

1ETA)*Y(24)/X-ETA*Y(2J*Y(20)/(S*X) 

221  70    RETURN 

222  ENO 

223  .         SUBROUTINE    CMINV< A,N,C, L,M) 

224  DIMENSION   A  (9)  ,L  (3)  ,!<{3  ) 

225  OOUBLE    PRECISION    A,C,B IGA ,HCLD,CABS        • 

226  0=1.0 

227  NK—N 

228  00  80  k=1,n 

229  NK=NK+N 

220  L1KJ--K 

221  M(K)*K 

222  KK-NK+K 
233  BIGA»A<KK) 
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234  00  20  J»K,N 

225  IZ*N*(J-1) 

226  00  20  I-K,N 

227  IJ»IZ+I 

226  10  IF(DABS<BIGA)-DABS(A<IJ])>  15,20,20 

239  15  BIGA*A(IJ) 

240  L(K)=I 

241  M»K)»J 

242  20  CCNTINUE 

243  J=L(K) 

244  IF(J-K)  35,35,25 

245  25  KI=K-N 

246  00  30  1*1, N 

247  KI=KI*N 

248  HCLD=-A(KI) 

249  JI=KI-K+J 

250  A<KI)=AUI) 

251  30  A(JI)=HCL0 

252  35  I=M(K) 

2*2  IF(I-K)  45,45,38 

254  38  JP=N*(I-l) 

255  00  40  J*1,N 

256  JK=NK+J 

257  J1=JP+J 

258  HCLC=-A(JK) 

259  A(JKJ=A(JI) 

260  40  A(JI}«hCL0 

2(1  45  IF(BIGA)  48,46,48 

262  46  0=0.0 

263  RETURN 

264  48  DC  55  1*1, N 

265  IF(I-K)  50,55,50 

266  50  IK=NK«-I 

267  AUK)=A(IK)/(-BIGA) 

268  55  CCNTINUE 

269  00  65  1=1, N 
27C  IK=NK*I 

271  HCLD-AUK) 

212  IJ=I-N 

273  00  65  J*1,N 

274  IJ=IJ*N 

215  IF(l-K)  60,65,60 

276  60  IF(J-K)  62,65,62 

277  62  KJ*IJ-I*K 

278  A<  IJ)=HCLO*A(KJJ-»-A(lJ) 

279  "    65  CONTINUE 

260  KJ=K-N 

261  DO  75  J*l,N 
282  KJ=KJ*N 

263  IF(J-K)  70,75,70 

264  70  A(KJ)=A(KJ)/BIGA 

265  75  CONTINUE 

266  D-C*BIGA 

267  A<KK)=«1.0/BIGA 

268  80  CCNTINUE 

289  K=N 

290  100  K=(K-i) 

291  IF(K)  150,150,105 

292  105  I=L(K) 

293  IF(I-K)    120,120,108 


Computer  program  for  annular  orthotropic 

convex  variable  thickness  plate  — 

Forward  shooting. 
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SJOB 

(>***»****  *******  4*  +*  +  «*.**»*<•  *************  *******  ******** 

C  INITIAL  VALUE  METHOD  -  FREE  VISPATICN  OF  AN  ANNULAR 

C  ORTHCTRCPIC  .CCNVEX  VARIABLE  THICKNESS  PLATE.  WITH 

C  BOUNDARY  CONDITIONS  AS  FIXEO  CN  THE  INS  ICE  ANO 

C  FREE  ON  THE  OUTS  ICE 

£*************»********  »*»»»**»»»**»  *************  ******** 

C** ************#************* ^♦v**************** ******** 

C  S*RATIO    CF    ELASTIC    CONSTANTS 

C  E=PCISSCNS    PATIO 

C  OL=UNIFORM    LOADING    INTENSITY 

C  A=AMPLITUOE 

C  R=RATIO    CF     INNER    TO    CUTER    RACIUS 

C  P=£IGENVALUE 

C  H=STEP    SIZE 

C  ETA=TbICKNESS  FUNCTION 

C  DETA=FIRST  CERIVATIVE  OF  ETA 

C  OOET=SECCND  DERIVATIVE  CF  ETA 

£*******  *******  *******  *•<**  ******  .<!**  ***<:*******  *********** 

1  IMPLICIT  REAL*3(A-H,0-Z>. INTEGER! I-N) 

2  DIMENSION  ETA141)  ,XX(4l) ,Y(3U),C(3U  )  ,TP  14  ,4  )  .  C  (8,  41  > 

3  DIMENSION  C(4).MU(4> ,LW(4) ,ER(4) 

4  OIMENSICN  R8SI45).CeS(45),RMS(45) .CKS145) 

5  112  FORMAT(5X.'AMP='  .022 . 14  ,3 X , • FREC= «  ,D22 . 14  ,3X, • FRER= • 

1022.14) 

6  113    F0RMAT(9X.'W'  ,  19X  . ' CW ' , 18X ,  'DCW,  17X,'0CCW> 

7  114    F0RMAT(4D22.14J 

8  115    F0RMATI//9X,  'F',  19X,  'OF'  J 

9  117    FORMAT! 1H    ) 

10  120    F03MAT(5X.'S=«,r-lC.3.5X.'E=',FlC.3,5X,'QL  =  '  ,F10.3> 

11  121    FCRMATI//9X.' STA',19X. 'PROF') 

12  12?  F0RMAT(6X.'R8S' .ISX.'CnS' ,18X. 'RMS'  ,18X, 'CMS' ) 

13  123  F0RMAT(5X, • ITER=«. 12) 

14  S=0.5 

15  H=l./40. 

16  Ll  =  41 

17  KK=9 

18  JK=LL  +  1--K< 

19  IK  =  l 

20  0L=0.0 

21  A=0.0 

22  6=1. /3. 

23  R=0.2C-0 

24  0A=0.25 

25  P=4.0**2 

26  ETA1=8.00 

27  ETA2=-58.0 

28  ETA3=0.06 

29  8E=0.5 

£********  *****-?************************************  ****** 

C  CONSTRUCT  INITIAL  VALUES 

C**#*  ******************************  ******  **************** 

30  5C0  IT=1 

31  DO  9  1=1.30 

32  9  Yt I)=0.00-0 
23  YCl)=l.OC-0 

34  Y(3)=ETAl 

35  YU)=ETA2 

36  Y(5)=ETA3 

37  Y!6)=E+ETA3/(R*SJ 
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60C 

43 

44 

623 

45 

46 

620 

47 

48 

4<5 

50 

51 

615 

52 

6?4 
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38  Y<9)=1.0 

39  YC16>=1. 

40  Y(23)=1.0 

41  Y(24)=E/(R>$) 

£**,*********  +  *>«**.******•*************♦******************* 

C  X*INCEPEN0ENT    VARIABLE 

C  INTEGRATION    BY    FCPWARD    SHCCTING 

C********** **************************** ****************** 

X=0.20-0 

00  622  1=1.30 

0(1  )=Q.CC-0 

OC  620  1=1.6 

C(  I,KK)=YI I) 

KJ=KK+i 

00  624  J=KJ.LL 

CALL  RKGU.H.Y.Q.P.A.S.E.QL.BE) 

DO  615  1=1.6 

D(L.J)=Y(L) 

CONTINUE 
C**4*»*¥*******  *.********«.***♦***♦*********>■***  *********** 
C      ER(I)=ERRO*  VECTOR  FOR  BOUNDARY  CCNCITICNS 

Q********** **************************************** ****** 

53  ER<  l)=CI  l.LD-l.U 

«4  ER<2)=E*C(2.LL)+S*D( 2.LL) 

55  ER(3)=S*0(4,LL)-S*0.5*0(3.LLJ-(l.*1.5*E)*0C2.LL) 

56  ER(4)=0(5.LL) 

C*** ********♦<******** **v ****************  **************** 

C      CONSTRUCT  ERRCR  NCR* 

C ******  ******  4  **********  ****•«****»<*******  **************** 

57  CO  26  1=1.4 

«S  DFR=OAeS(ER(I) J 

59  1F(DER.GT.0.1D-G5)  GC  TO  23 

60  26    CONTINUE 

61  GO    TC    900 

62  28    CGNTINUE 

C**************** ****************** ********************** 

C  NEVTONS    HETHOO    (EPRCR    NCRM=(B)Y=0 

C  TP(I,J)=    THE    JACCSIAN   OF    THE    INITIAL    VALUES 

C  C(I)=CCRRECTICN    VECTCR 

C****** ******************************  ******************** 

63  TP(i.l)=Y(7) 

64  TP(2.1)=E*Y(8)*S*Y(9) 

65  TP(3,l)=-(l.*1.5*£]*Y(a)-S*0.5*Y{9)+S*Y(10) 

66  TP(4.l)=Y(ll) 

67  TP(1,2)=Y<13) 

68  TP(2.2)=E*Y(14)+S*Y( 15) 

t<i  TP«3.2)=-(1.-H.5*E)*Y(14)-S*0.5*Y(15)+S*Y(16) 

70  TP(4,2)=Y( 17) 

71  TP(1.3)=Y(19) 

72  .     TP(2.2)*E*Y(20)*S*Y(2U 

73  TP<3.3)=-(l.+l.5*E)*Y(20)-S*0.5*Y(21)*S*YC22) 

74  TP(4.3)=Y(23) 

75  TP( 1.4)=Y(25) 

76  TP<2.42=E*Y(26)+S*Y{27) 

77  TP(3.4)=-(i.*1.5*E>*YI26)-S*0.5«-Yl27)*S*Y(23) 

78  TP(4.4)=Y(2SJ 

79  DET=O.0C-O 

€C  CALL    Of  INV{TP.4.CET.LW,HW> 

81  DO    75    1=1.4 

82  C  I  I  1=0.0 
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£3  OC    75   J=1.4 

84  75   CU)=C{  !)-TP(  I.J)*ER(J5 

C*»  *»**»#*****««  •****#***»»***#*♦*»»*  *»%*♦»************** 

C  CCRRECTcC    VALUES 

C** ********** ****** ******************************** ****** 

£5  OQ    76    5=1.6 

£6  76    Y(I)=C( r,KK) 

£7  Y(3J*Y(3)*CU) 

£6  Yt4)=Y(4)+C(2) 

e?  Y(5)=Y(5)+C<31 

90  P=P+C(41 

?l  YC6)=E«Y(5)/(R*S) 

"52  DO   77    1=7,30 

S3  77    Y(I)=O.0D-O 

94  Y<9)=1.0 

95  Y(16)=l.0 

96  Y(23)=1.C 

<57  Y(24)=£/(R*S) 

98  IT=IT*l 

$<}  IFl  IT.GT.IC)    GC    TO    550 

ICO  GO   TO    600 

C* *♦*♦*****+********+*****************+ ****************** 

C  FINAL    RESULTS 

C************ ******************************** ************ 

101  9G0  SRA=OSCRTtA) 

IC2  SP=0SCPT(P) 

IC3  DO  795  J=KK,LL,4 

IC4  0J=J-1 

105  XX(J)=CJ*H 

1C6  ETA(J)  =  L-8E*(XXU)+*  10.5)) 

107  IF(XXU).GT.O.O)  GO  TO  905 

1C8  RBSU)=SRA*C13.J)/2.*(1.-EJ 

109  CBS( J)=RBS( J) 

HO  RMS(J)  =  A*C<6tJ> 

HI  CMS(  J)=R^S(J) 

112  GO  TO  795 

113  905  R8S(J)=SPA*ETA{J)*(S*0(3,J)+E*0(2.J)/XX( J))/2.*(S-E**2) 

U4  C8S(  J)=SR.VrETA{  J  )*  (0  (2  .  J  )  /  XX  {  J  )  +£*C(3,J)  )/2  .*  (  S-E*  *2  ) 

115  RMS(J)  =  A*0<5,J)/ETAU)*XX(  J) 

116  CKS(J)=A*C(6.J)/ETAU) 

117  795  CCNTINLE 

C*********** ***********  ********************************** 

C  FCR  FRECLENCY  RATIO 

<>♦*+** ****** **************************************  ****** 

118  IFtA.GT.O.O)  GO  TC  906 

119  SPC=SP 

120  9C6  SPR=SP/SPO 

121  WRITE(6,U7) 

122  WRITE16.120)  S.E.CL 
12*3  WRITE(6.117) 

124  WRITEI6.112)  SRA.SP.SPR 

125  WRITE(6.117) 

126  WRI7EJ6.U3) 

127  00    901    J=KK,LL,4 

128  901  WRITEI6.114)  {0 ( I  , J )  .1  =  1 14 ) 

129  WRITE(6,ll5) 

130  OC  902  J=KK.LL.4 

131  9C2  WRITEC6.ll*)  (D( I  , J)  »L=5,6) 

132  WRITE(6.117) 

133  WRITE16.122) 
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134  00  903  J=KK,LL.4 

135  9C3  WR1TEI6.1U)  PBS  {  J  )  .C3S  I  J)  .  RPS  I  J)  t  Cf  S  ( J  1 

136  WRITEU.IW) 

137  WRITE16.123)     IT 

138  WRITEI6.121) 

129  00    924    J=KK,LL.4 

140  924    WRITEC6.U4)    XXU).ETA(J) 

141  WRITEC6.117) 

<>¥****+* ******************************** **************** 

C  PERTUReATlCN    CF    AMPLITUDE 

C**+*********»»?***+v*  *****•****************************'''* 

142  A=A*CA 

143  IK=IK+l 

144  IFUK.GT.26)  GO  TO  550 

145  ETAl=C(3.KK) 

146  ETA2=C(4.KK) 

147  ETA3=C(5.KK) 

148  "  P*(SP-0.3)**2 

149  GO  TO  500 

150  550  STOP 

151  END 

152  SUBROUTINE    RKG( X .H. Y.Q. P. AP .S tE  .OL  ,  BE) 

153  IMPLICIT    REAL*8(A-H,0-Z), INTEGER!  I-N) 

154  Dlf-ENSICN   Y(30) .C(30) »DY(30).A(2) 

155  At 1)=0.292893218E135 

156  A12)=1.7C7106781L£65 

157  H2=0.5*H 

158  CALL    CERIV(X.h.Y,DY.P.AP.S.E,QL.BE) 

159  00    13    1=1.30 

160  R=H2*0Y(  1 1  —  0 ( I ) 

161  YtI)=Y(I)*R 

162  13   0<l)=OUJ+3.0*R-H2*OY(I) 

163  X=X*H2 

164  00    60    J=  L  .2 

165  CALL    DERIVIX.H.Y.OY.P.AP.S.E.CL.BE) 

166  00    20    1=1.30 

U7  R  =  AtJ)*(H*OY(I)-C(in 

168  Yd)»YtI)+R 

169  20    0(I>=C(  1)4-3. 0*R-A(J)*H*CY(1) 
17C  60   CONTINUE 

171  X=X+H2 

172  CALL    CERIV(X,H,Y.0Y.P.AP.S.E,CL.8E) 

173  DO    26    1=1.30 

174  R=(H*CY(  I)-2.0*0(  I)  1/6.0 

175  Y(I)=Y(IJ*R 

176  26    0(  I)=C(  IH-3.0*R-H2*0Y(I) 

177  RETURN 

178  ENO 

179  SUBROUTINE    CER  IV(  X.H. Y.OY , P. AP . S .E  ,CL  ,3E) 

180  IMPLICIT    REAL*8(A-H.0-Z)  ,  INTEGERd-N) 
161  DIPENSICN    Y(30).CY(30) 

182  ETA=l.-Ec*(X*v(0.5) ) 

183  0ETA=-C.5<«8E*(X**(-0.5)) 

184  D0ET=G.25*BE*(X**(-1.5)) 

185  00  10  1=1.3 
166  10  0Y( I)=Y(H-l) 
1£7  DY(5)=Y(6) 
188  00  12  1=7.9 
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169  12    0Y(  I)-Y(H-i) 

190  0YC11)»Y(121 

191  DO    15    1=13.15 

192  15   DY(  I  )=Y(I+i) 

193  OY(l7)=Y(13) 
ISA  DO    16    1=19,21 

195  16    DYUJ=Y(  IU) 

196  0Y(23)=Y(24) 
1^7  DO    20    1=25.27 

198  20   DY(l)=Y( 1+1) 

199  0YJ29)=Y(30) 

2C0  IF(X.C-E.C.lC-02)    GC    TO    17 

201  0YI4l«J.*P*V(l>/8.META"*2>+27.'*tl.-E**2)*AP*Y{3>*Y(6) 

l/4.*(ETA*-3) 
2C2  DY14)=CY(4)-9.M 1  .  +  £ )*DOET*Y( 3 > /8 ,*ETA 

203  1F(AP)     18.18.19 

2C4  19    DY(4}»OY(4l+3.*(l.-E«*2)*OL/8.*(ETA**3l*0SGRT(AF) 

205  18    DY(4>=0Y{4)  ,     ._,_ 

206  0Y( 10>  =  3.^P*Y(7)/S.*(ETA**2H-27.*C1.-E**2)*AP*(Y(9)  + 

1Y(6)>Y(3)^Y(12) }/4.*(ETA**3) 

207  DY<l0)=CY(10)-9.*(l.*EKDCcT*Y(9)/8.*ETA 

208  DY(  16>=2.*P*Y(13>/8.META**21+27.*(l.-E-*2)*APi'<Y(i5)* 

1Y(6)*Y( l£)*YU) )/4.*lETA**3) 
^C9  0Y(16J=0Y(16)-9.Ml.fE)*0DET*Y(l5)/8.*ETA 

2io  0Y(22  1«2.-fP*Y(19J/8.*(ETA**2J*3.*YCl)/a.*CETA**31+27. 

1*11.-E**2)*AP*(Y(21J*Y(6)+Y(3)*Y124>)/4.*(ETA**3) 

211  0Y( 22 ) =0 Y I  22 ) -9. *< I. * E) *00E T* Y( 21)/ 8. *ET A 

212  DY(6)=0.0 

213  DYti2)=0.0 

214  0Y(18)=0.0 

215  0Y(24)=0.0 
?lfc  0Y(30)=0.0 

217  GO    TO    70 

218  17    80=19. *(S-IE**2)  >)/tS*X*(ETA**3)) 

219  Bl  =  (S-E**2J/(  (l.-E**21"(S<META**2))) 

220  32=U6.*CETA)/£TA)t(2./XJ 
B3  =  <-1./(S"U**2)))*(3.*0ETA*(2.*S  +  E) )/ IS*X*ET A ) *< 3 .* 

10DET>/ETA+t6.*(OETA**2)>/<ETA**2) 
d4=l./(S*(X**3)  )-(3.*0ETA)/(S*(  X**2  )*ETA  )  +  (  3.  *E*0DtT  )  / 


221 
2  22 
223 


l(ETA*S*X)*-G.*E<-(CfcTA+*2)/{S*X*(ETA:«  +2)  ) 

0Y(4)=SG»AP-<YC3)*Y(5)+Y(2)*Y(6))+ai*P*Yll)-B4*Y(2) 
l-B3*Y(3)-B2*Y(A) 

224  IF(AP)     1<_U. 100,101 

225  101    DYl41»CY(4)  +  <S-lE**2n*QL/(S*(ETA**3>*0SQRTlAP») 

226  100    DY(4)=0Y(4J  .«-..,-*. 

227  OYI6)»U.-E»X*CETA/ETA)*Y(5)/CS*(X**2)l-(l--X*0ETA/ETA 

i)*Y(6)/X-{ETA*(Yl2)**2)  )/  <2.*X*S> 
DY(  lo)«B0*AP*(Y(6)*1f(aj*Y(5J*Y(S)+Y(2}*Y(12)*Y(3)*Y(ll 

1)  J481*P*Y( 7)-82*Y(  10  5-B3*Y(9  3-B4*Y( 8) 
OY(l2l»Cl.-E*X-*0ETA/ETA)*mi)/lSMX**2>»-U.-X*OETA/ 

1ETA)»Y(  12)/X-ETA*Y(2)*Y(8)/(S*X) 

0Y(loJ=eC*AP*(Y(5)*Y( i5)+Y(3>*YU7)*YC  14  )*Y  (  6) +Y  (  2)  * 
1YC 18) )+ei*P*Y<  13)-02*Y(L6)-83*Y(15)-84*Y{14) 

DY(18)=(1.-E«X*0ETA/CTA)*Y< 17 )/ (S*( X**2) )-( 1. -X*0ETA/ 
IETA)*Y{18)/X-ETA'»Y(2)*Y(14)/(S*X) 

DY<22>=BG»AP"(Y(21)-Y(5)*Y<3)*Y(23)«-Y(20)*Y(6) *Y(2)* 
1Y(24)  )  +  Pl*P*Y< 19)+B1*Y(  l)-B2*Y(22)-83*Y(2i)-3  4*Y(20) 

DY124)=( l.-E*X*0£TA/ETA)*Y(23 )/ tSMX**2) )-( l.-X*D£T A/ 
1ETA)*Y(24)/X-ETA*Y< 2  )  *Y ( 20 ) / ( S* X) 

DYl28)  =  BG*AF*(Y(27)*Y(5)*Y(3}*Y(29)*Y(6)+Y(26H-Yt2>* 


228 
229 
230 
231 
222 
233 
234 
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lY!30))*8l  +  P*Y(2  5)-a2*Y(28)-B2*Y<27)-B4*Y(26}*-Bl*Y(l) 
235  DY(30)=(1.-E*X«0ETA/ETA)*Y129>/{S*(X**2)) 

l-d.-X*  CETA/ETA)  ?Y(3C)/X-ETA*Y(26)*Y(2)/(S*X) 


236 

70 

RETURN 

227 

END 

228 

SUBROUTINE  OMI  NV  ( A  ,N  .Of  It  »*) 

239 

DIMENSION  A(16) ,L(4).M(4) 

240 

0CU8LE  PRECISION  t .0 , 8IGA ,HClC, CAES 

241 

0-1.0 

242 

\K»-N 

243 

00  80  K=1,N 

244 

NKsNK+N 

245 

L(K)=K 

246 

M(K)*K 

247 

KK=NX*K 

248 

BIGA=A(KK) 

249 

00  20  J  =  K,N 

2  50 

IZ=N*(J-1) 

251 

00  20  I=K.N 

252 

IJ=IZ*I 

253 

10 

iF(CAes(eiGA)-0Aes(A{U)  j  )  15,2c, 20 

2  54 

15 

BICA=A( IJ) 

255 

UK)=»I 

256 

MCK)=J 

257 

20 

CONTINUE 

258 

J=l(K) 

259 

IF(J-K)  25,35,25 

260 

25 

KI=K-K 

261 

00  30  I=1,N 

262 

KI=KI+N 

263 

HOLO=-AIKI) 

264 

JI=KI-K+J 

265 

A(KI J=A( JIJ 

266 

30 

AUI)=HCLO 

267 

35 

I=M(K) 

266 

IF(I-K)  45,45,38 

269 

38 

JP=N*(  I-l) 

270 

00  40  J=l,N 

271 

JX=NK*J 

2  72 

JI=JP+J 

273 

HCL0=-A( JK) 

274 

AUK)  =  A!JI  ) 

275 

40 

AUI)=HQLD 

276 

45 

1FIBIGA)  43,46,48 

277 

46 

0  =  0.0 

278 

RETURN 

279 

48 

00  55  1=1, N 

2  80 

IF(I-K)  50,55,50 

2-81 

50 

IK=NK*I 

262 

A(IK)=M IK)/(-8IGA) 

2E3 

55 

CGNTINUE 

2  84 

DC  65  1=1. N 

265 

IK=NK+I 

2  86 

HOLO=A(IK) 

2e7 

IJ=I-N 

288 

00  65  J=1.N 

289 

1J=IJ*N 

290 

IF(I-K)  60,65,60 

291 

60 

IFU-KJ  62,65,62 

63 


. 


292 

62  KJ*IJ-I*K 

2  93 

A<  U)*hCL0*A(KJ)+A(IJ) 

294 

65  CONTINUE 

295 

KJ=K-N 

296 

DO  75  JM.N 

297 

KJ=KJ+N 

29a 

IF(J-K)  70.75,70 

299 

70  A{KJ)=A(KJ)/BIGA 

3C0 

75  CGNTINLfc 

301 

0=C*BIGA 

3C2 

A{KK)=1.0/BIGA 

3C3 

80  CCNTIME 

304 

K=N 

3C5 

100  K=(K-1) 

306 

1FU)  150,150.105 

307 

105  I=L(K) 

3C8 

IFCI-K)  120,120.108 

309 

108  J0=NMK-1) 

310 

JP=N*U-1) 

311 

DO  110  J=l.N 

312 

JK=JO*J 

313 

HOLD-A(JK) 

31* 

JI*JR«-J 

315 

AUK)=-A(JI  ) 

316 

110  A{ JI >=KID 

317 

12U  J=MK) 

318 

IMJ-KJ  1U0.100.125 

319 

125  KI=K-N 

320 

DO  130  1=1, N 

321 

KI=KI+N 

322 

HCLD=A(KI) 

323 

JI=KI-K*J 

324 

A(K1)=-M  JI) 

325 

130  A( JI)=HCLD 

326 

GC  TO  ICO 

327 

150  RETURN 

326 

FND 
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ABSTRACT 

The  problem  of  finite  amplitude,  axisymmetric  free  vibration  of 
variable  thickness  orthotropic  annular  plates  is  formulated  in  terms  of 
the  Von  Karman's  dynamic  equations.  A  kantorovich  averaging  technique 
is  applied  to  convert  the  nonlinear  boundary  value  problem  into  a  cor- 
responding eigenvalue  problem  by  elimination  of  the  time  variable.  A 
numerical  study  is  proposed  by  introducing  the  related  Initial  value 
problem.   By  making  successive  corrections  and  perturbations  of  the 
parameters  in  a  numerical  solution  to  the  initial  value  problem,  approxi- 
mate solutions  to  the  boundary  value  problem  are  ohtained.   The  cases 
investigated  are  free  outside  and  fixed  Inside,  parabolic  and  convex 
variable  thickness  orthotropic  annular  plates. 

The  hard  spring  behavior  is  evident,  and  it  Is  found  that  the  mode 
shape,  bending  stresses  and  membrane  stresses  are  nonlinear  functions  of 
the  amplitude  of  vibration. 


